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Friction is one of the main sources of dissipation at liquid water/solid interfaces. Despite recent 
progress, a detailed understanding of water/solid friction in connection with the structure and 
energetics of the solid surface is lacking. Here we show for the first time that ab initio molecular 
dynamics can be used to unravel the connection between the structure of nanoscale water and 
friction for liquid water in contact with graphene and with hexagonal boron nitride. We hnd that 
whilst the interface presents a very similar structure between the two sheets, the friction coefficient 
on boron nitride is « 3 times larger than that on graphene. This comes about because of the greater 
corrugation of the energy landscape on boron nitride arising from specific electronic structure effects. 

We discuss how a subtle dependence of the friction on the atomistic details of a surface, that is not 
related to its wetting properties, may have a significant impact on the transport of water at the 
nanoscale, with implications for the development of membranes for desalination and for osmotic 
power harvesting. 


Nanofluidics is an exciting field that offers alterna¬ 
tive and sustainable solutions to problems relating to en¬ 
ergy conversion, water filtration and desalination m- 
Miniaturization towards nanofluidic devices inevitably 
leads to an enhanced influence of surface and interface 
properties as opposed to those of the bulk. Friction is the 
most important interface property that limits fluid trans¬ 
port at the nanoscale, and its understanding is therefore 
crucial for the design of more efficient membranes, nan¬ 
otubes and pores that exhibit low liquid/solid friction. 
The behavior of liquid flow at scales on the order of a 
few tens of nanometres departs from continuum fluid dy¬ 
namics and desirable transport properties emerge at such 
small scales m- For instance, carbon nanotubes have a 
very high water permeability as compared to the predic¬ 
tion of macroscopic fluid dynamics (2]. Further, a van¬ 
ishing friction has been found, giving rise to superlubric 
behavior of water chains inside tubes of sub-nanometre 
radii [TT] . 

Besides carbon, boron nitride (BN) nanostructures 
have recently been explored for the development of 
nanofluidic devices for fast water transport and efficient 
power generation [DII1II3]. Recent interest has been 
fueled by the demonstration that salinity concentration 
gradients across BN nanotube membranes can leed to the 
generation of very large electric currents [1]. It has also 
very recently been shown that there is a very large inter¬ 


layer friction between in multiwalled BN nanotubes [14j . 
as opposed to the superlubric behavior of the (homopo- 
lar) carbon nanotubes [T^. This suggests that the fric¬ 
tional properties of BN and C nanostructures might be 
quite different. However, to the best of our knowledge 
there has been no attempt to measure or compute the 
friction of water at the interface with BN sheets or nan¬ 
otubes. Given that ab initio results have shown very 
similar contact angles of water droplets on graphene and 
BN sheets it remains to be seen if transport proper¬ 
ties on these two systems are also similar. 

The rise of the atomic force microscope and the sur¬ 
face force apparatus has advanced our understanding of 
nanoscale liquid/solid friction substantially [17]. Yet, it 
remains extremely hard to relate friction and dynamics 
to structure and wetting of solid surfaces from experi¬ 
ment. This is especially true for the case of graphene 
where there is even controversy over the water contact 
angle (see e.g. Refs. [T8l[T9]). Molecular simulations, and 
especially force field molecular dynamics, have proved ex¬ 
tremely useful in investigating structural and dynamical 
properties of confined liquids elucidating molecular level 
information that is challenging for experiments to ob¬ 
tain. Because of the dependence of simulation results 
on the parametrization of force fields, it is interesting to 
investigate the dynamics of interfacial water using elec¬ 
tronic structure methods so that two different materi- 
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als such as graphene and BN can be compared on an 
equal footing. Accordingly, ab initio molecular dynamics 
(AIMD) offers an interesting alternative that has been 
widely used to study complex liquid/solid interfaces (see 
e.g. Refs. [121 HOI HI] ) • Although transport properties 
in liquids have been computed before using AIMD (see 
e.g. Refs. we show for the first time that it is 

possible to compute a converged friction coefficient from 
AIMD. 

In this work we compare the structure and dynamics 
of water in contact with BN and carbon nanostructures 
from AIMD. Specifically, we study liquid water in con¬ 
tact with graphene and hexagonal BN sheets, which is 
relevant to understanding flow at membranes based on 
these materials and also inside large nanotubes [25|. We 
find a striking similarity between the structure of water 
in contact with the two sheets. Nevertheless, there is a 
three-fold increase in the friction coefficient on BN be¬ 
cause of a more corrugated free energy surface on BN 
compared to graphene. This work illustrates the com¬ 
plexity of nanoscale friction where subtle electronic ef¬ 
fects, with no detectable consequences on the structure 
of the liquid, may however have a dramatic impact on 
water transport at the nanoscale. 

We performed a series of extensive AIMD simula¬ 
tions of a thin liquid water film on graphene and on a 
single layer of hexagonal BN. Full details of the com¬ 
putational set-up can be found in the Supporting In¬ 
formation (SI). However, the key features are that we 
used the CP2K [5^1 code to perform AIMD simula¬ 
tions in the Born-Oppenheimer approximation with the 
electronic structure computed at the density functional 
theory (DFT) level using the optB88-vdW exchange- 
correlation functional d?] HHj- The optB88-vdW den¬ 
sity functional has been shown to predict accurate in¬ 
terlayer binding and interlayer distances in graphite and 
bulk hexagonal BN [29] , to give a good description of the 
structure of bulk liquid water [3D] and also to correctly 
describe the relative stability of water-ice structures on 
metals compared to the bulk ice lattice energy [3T]. The 
water/graphene and water/BN sheets have been mod¬ 
elled using 6 X 10 orthorhombic cells about 25 x 25 
wide, with « 20 A thick liquid water films. There is 
a vacuum gap between the liquid-vacuum interface and 
the next periodic image of « 15 A. Each film contains 
400 water molecules and in total each system consists 
of 1440 atoms. Upon these systems we performed 40 ps 
long AIMD simulations with the last 35 ps of each tra¬ 
jectory used for analysis. As an illustration snapshots 
from the AIMD simulations for the water/graphene and 
water/BN interfaces are shown in Fig. [^c) and (d), re¬ 
spectively. The AIMD simulations are in the canonical 
ensemble close to room temperature (330 K). We also 
performed a number of additional AIMD and force field 
MD simulations to explore the sensitivity of our results to 
the use of a different ensemble and to issues such as finite 


size effects, different initial conditions and time scales. 
Overall, the results of these tests, which can be found in 
the SI, agree with the results presented herein. 



FIG. 1. Structure of the liquid water film on graphene (GRA) 
and a single sheet of hexagonal BN. Average density profile 
(p) (a) and number of hydrogen bonds (H-bonds/H20) (b) as 
a function of the height from the sheet. In (b) the geometric 
criterion from Luzar and Ghandler was nsed to dehne a H- 
bond |32| . Snapshots of the liquid him on GRA (c) and on 
BN (d). In (c) and (d) O and H atoms are coloured in red and 
white, while G atoms are in light blue, and B and N atoms 
are in pink and purple, respectively. The liquid him structure 
is essentially the same on both sheets, only on BN the him 
is slightly thinner because the BN unit cell is ~ 3.5% larger 
than that of graphene. 

We begin our analysis by illustrating in Fig. [^a) the 
planar average density profile (p) as a function of the 
height from the sheets. The graphene and BN sheets 
exhibit oscillations in the heights of the atoms within 
them of up to 1.2 A. To account for this, we compute 
the height of an atom in the liquid water overlayer as 
the height difference to the closest atom in the sheet. By 
computing the density profile in this way we fully account 
for the layering of the liquid film, which would otherwise 
be partially hidden behind the oscillations of the sheets. 
The two density profiles on graphene and BN overlap for 
most of the film height, which is the first signature of the 
apparent similarity of the structure of liquid water on the 
two sheets. Perturbation from the bulk liquid induced by 
the surface is most significant within the first 10 A of the 
surface. This is consistent with previous reports on liquid 
water/solid interfaces (see e.g. Refs. [DO] HI] |33| [51]1. 
which extends in this case to about 8 A from the sheets. 
Within this region there are two evident peaks for each 
of the two density profiles. The first at a height of about 
3.0 A hits a density maximum of « 3.7 g/cm^. After this 
first peak there is a depletion of water with a minimum 
at about 4.5 A with a density of 0.4 g/cm^. We define the 
contact layer as the region delimited by this minimum in 
the density profile. At « 6.0 A the second peak appears 
with a density of « 2 g/cm^. Further away from the 
sheets, density oscillations are gradually suppressed and 
the density of liquid water is recovered with an average 
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value ofl.l2±0.16 g/cm^. Finally, the density decays 
as it is characteristic of the liquid-vapor interface [35]. 
The decay in the density on BN starts off at about 17 A, 
slightly before than on graphene simply because the BN 
unit cell is « 3.5% larger. 

As further confirmation of the striking similarity be¬ 
tween the structure of liquid water on graphene and 
BN, we compare the profile of the average number of Id- 
bonds per molecule, (H-bonds/H20) in Fig. Bb). Also 
in this case the two curves essentially lie on top of each 
other for all of the film height apart from the liquid- 
vapour interface region. The steep increase in the (H- 
bonds/H 20 ) starting at about 2 A shows that already 
in the contact layer water molecules engage in a large 
number of H-bonds of « 3 /H 2 O. Shortly after this initial 
increase, there is a drop to « 2.5 H-bonds/H20 corre¬ 
sponding to the depletion region around 4.5 A from the 
sheet. The pronounced fluctuations in the number of wa¬ 
ter molecules in this region partly penalizes H-bonding 
with neighboring waters. After the depletion region, at 
a height of 5 A, the number the (H-bonds/H20) rises 
to the bulk value of « 3.5. It remains constant until 
the liquid-vapour region is approached, at a height of 
17 A, where there is a rapid decay to zero in the (H- 
bonds/H 20 ) within 3 to 4 A. Finally, other structural 
characteristics of the two systems exhibit striking simi¬ 
larities, such as the orientations of the water molecules 
within the films (see Fig. S2 in the SI). 

Having compared the structure of liquid water on the 
two sheets we now turn our discussion to investigate its 
dynamics. Specifically, we focus on the friction coefficient 
A, defined as the ratio between the friction force parallel 
to the sheet Fp per unit area (A) and VsUp, the velocity 
jump at the interface m, namely A = Fp/{AvsUp). In 
the framework of linear response theory, A can be ob¬ 
tained from the equilibrium fluctuations of the friction 
force, using a Green-Kubo relation |31I37|: 

A = lim XGK{t), (1) 

t—¥00 

with 

where fee is the Boltzmann constant, T the temperature 
and the factor of 1/2 comes from the averaging over the 
two spatial dimensions parallel to the sheets. We show in 
Fig. AGif(t) for the case of water on graphene and on 
BN, which for sufficiently long time intervals (here ap¬ 
proximately > 0.3 ps) plateaues to the value of A. The 
key result is that the friction coefficients on the two sheets 
differ significantly, with the friction on BN being approx¬ 
imately 3 times larger than on graphene. Specifically, 
while A = (9.6 ± 2.0) x 10^Nsm“^ for graphene, we ob¬ 
tain a value on BN of A = (30.0 ± 5.4) x lO'^Nsm”^. A 
number of tests using both AIMD and force field MD, 


ensured that our results on the friction coefficients are 
converged and that the increase of about three times in 
the friction on BN compared to graphene is consistently 
observed, regardless of the type of ensemble used, time 
scales, or different system sizes. The results of these tests 
are presented in Fig. S3 in the SI. 

Whilst liquid/solid friction is the relevant microscopic 
property that quantifies the dynamics of a fluid at the 
nanoscale, a length scale characteristic of the flow is of¬ 
ten measured experimentally. This is the so called slip 
length b, defined as the distance relative to the surface 
where the linear extrapolation of the tangential flow ve¬ 
locity vanishes. We can relate the slip length b to the 
friction coefficient via the shear viscosity of bulk liquid 
water 77: b = ij/X |10j . Using the experimental bulk wa¬ 
ter viscosity rj = 10“^ Pas, the corresponding slip lengths 
for graphene and BN are 10.4 ± 2.2 and 3.3 ± 0.6 nm, re¬ 
spectively. Overall water slippage on graphene and BN is 
characteristic of hydrophobic surfaces with a low friction 
coefficient, while on hydrophilic surfaces such as mica, sil¬ 
icon or graphene oxide slippage is significantly inhibited 
with sub-nm slip lengths [351135] • 



time [ps] 


FIG. 2. Comparison between the Green-Kubo estimate of the 
friction coefficient of liquid water on graphene (GRA) and on 
BN. The shaded areas represent the uncertainties obtained 
by performing a block average. The friction coefficient A is 
given by the plateau value at long times. There is an evident 
increase in the friction coefficient on BN. 

To rationalize the difference in the two friction co¬ 
efficients in terms of structural and energetic contri¬ 
butions, we compute the free energy profile of water 
within the contact layer AG{x, y), defined as AG{x, y) = 
—kBT\nPo{x,y). Here, Po{x,y) is the spatial probabil¬ 
ity distribution function of the O atom of water within 
the contact layer at a point (x, y) projected onto the 
primitive graphene and BN unit cells. A similar approach 
based on the 2-D spatial probability distribution function 
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has been used previously to understand water slippage on 
model MgO surfaces [40]. Fig. [^illustrates the free en¬ 
ergy profiles resulting from this analysis. We notice the 
very small energy scale within kgT at room temperature 
indicating a mobile contact layer. Although, the aver¬ 
age liquid structures in the two systems are very similar, 
the free energy profiles for the motion of water in the 
contact layer exhibit some clear differences. First, the 
free energy minimum on graphene is around the hollow 
site with the maximum around the top-C site, in agree¬ 
ment with previous work [20j . The minimum on BN is 
for the oxygens to sit around the B-site, as well as on the 
hollow site, while the maximum is on the N site. Sec¬ 
ond, and more importantly, the BN free energy profile 
(Fig. ib)) is more corrugated than that of graphene: 
the maximum corrugation of the free energy is only 13 
meV on graphene, but it is 21 meV on BN. Although we 
are discussing very small energies, this « 60% increase in 
the corrugation is observed consistently in all our various 
AIMD simulations of these systems. Indeed, the corru¬ 
gation is already converged if the free energy profile is 
extracted from half the AIMD trajectory, as opposed to 
the full trajectory (see Fig. S5). As we now discuss, this 
increased corrugation is the main reason for the observed 
increase in the friction coefficient on BN. 

It is interesting that such small energy differences, fess 
than 10 meV can contribute to a three fold increase in the 
friction coefficient and it is worth looking into this more 
in detail. It is known that the leading term in the friction 
depends quadratically on the corrugation of the potential 
energy surface felt by the water in the contact layer AF, 
such that A ~ AV^ (see e.g. mm)- We can approximate 
the corrugation of the potential energy AV with the cor¬ 
rugation in AG as obtained from the maximum value in 
free energy profiles in Fig. i such that A ~ AG^. To test 
the validity of this scaling relation we computed the ra¬ 
tio between the friction of water on BN and on graphene, 
Abat/Agka and compared it with the ratio between the 
square of the free energy extracted from Fig. [^a) and 
(b). From the ratio between the friction coefficients we 
obtain Xbn/^gra = 3.1 ± 0.8. The corrugation of the 
free energy on BN is 21 meV, while it is 13 meV on 
graphene, such that the ratio {AGbn/AGcra)^ ~ 2.6, 
is within the error of the ratio between the two values of 
the friction. 

The free energy corrugation depends on the atomic 
and electronic structure of the surface and on the H- 
bonds that form at a specific interface. It will differ 
but be related to the monomer potential energy surface. 
By performing an extensive series of calculations (about 
4000 wave function optimizations per sheet) we mapped 
the potential energy surface of a water monomer on the 
sheets (see Fig. SI), to explore if the potential of a wa¬ 
ter monomer on BN is more corrugated than that on 
graphene. This confirms that the free energy of the wa¬ 
ter contact layer on BN is more corrugated than that on 



FIG. 3. Free energy profile of water within the contact layer of 
the liqnid projected onto the graphene (a) and BN (b) primi¬ 
tive nnit cells. Although the free energy profiles are relatively 
smooth, a larger corrugation is present on BN and some dif¬ 
ferences are also observed in the topography. Transparent C 
and BN atoms are superimposed on the contour plots. Most 
stable (c) and less stable (d) configuration for a single wa¬ 
ter monomer adsorbed on graphene. Most stable (e) and less 
stable (f) configuration for a single water monomer adsorbed 
on BN. Only a small part of the unit cells used for the cal¬ 
culation of the monomer adsorption calculations is shown in 
figures (c) to (f). 


graphene because the potential of the water monomer is 
also more corrugated. 

From the analysis of the different structures of a water 
monomer on the two sheets, we have identified the main 
sites that give rise to the most and least stable configu¬ 
rations of water on graphene and BN. The most stable 
structure on graphene is for water on the center of the 
hexagon with its dipole pointing down towards the sheet 
(Fig. i c)). On BN water preferentially adsorbs in the 
center of the hexagon with the dipole moment parallel 
to the sheet and with one of the 0-H bonds pointing to¬ 
wards the N atom (Fig. j^e)). Upon considering water 
monomer adsorption at less favorable sites (see Figs.j^d) 
and (f)) we find that indeed the potential energy surface 
of the water monomer on BN is more corrugated than 
that on graphene (see Table SI). To understand the rea¬ 
son for the different corrugation between the two sheets, 
we performed a decomposition of the interaction ener¬ 
gies of several water structures on the sheets at the more 
stable sites (on the center of the hexagon of graphene 
and BN) and at less stable sites (on the top-C and top- 
N site), as shown in Table SI. Although there is a cer¬ 
tain level of arbitrariness to most decomposition schemes. 
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we have used an established approach to gain qualitative 
and semi-quantitative insight into our adsorption systems 
(see e.g. Refs. [29l[3l]). This analysis reveals that van 
der Waals dispersion forces and local-correlation effects 
do not contribute to the corrugation of the graphene and 
BN sheets. On the other hand, the interactions arising 
from exchange, electrostatic and electron kinetic energy 
terms are overall more pronounced on BN and provide 
the reason for the increased corrugation, and thus for the 
3-fold increase in the friction. Finally, we note that even 
though very small energy differences are being reported 
for the corrugation of the free energy profile, our AIMD 
simulations are long enough that the free energy profiles 
are converged even if they are extracted over half the tra¬ 
jectory (see Fig. S5). Also, these differences are robust 
with respect to the accuracy of the underlying potential 
energy surface (see Fig. SI) and they are consistent with 
the monomer interaction energies computed in Table SI 
using different codes and exchange-correlation function¬ 
als. 

Since this is the first AIMD study of water slippage at 
either of the two sheets examined, we wish to compare 
our results with the existing literature. In the absence of 
experiments on graphene or BN we compare our results 
to those obtained from atomic force microscopy measure¬ 
ments of water droplets on atomically smooth highly or¬ 
dered pyrolitic graphite, which have given slip lengths of 
between 8 and 12 nm [351112) in good agreement with 
our value for the slip length of 10.4 ± 2.2 nm. Although 
on graphene water transport may differ from that on 
graphite due to e.g. screening effects that may influence 
the binding on the two systems |43j . we do not expect 
the friction to change between graphene and graphite by 
more than a few percent. 

Compared to previous simulation studies, we notice 
that there is some uncertainty over the slippage of water 
on graphene obtained from force field MD, with values for 
the slip length between 1 and 80 nm [33] , and our AIMD 
result falls in this somewhat large range. The value of 
the friction on graphene obtained from our force field 
MD is « 3 X lO^Nsm”^ (see Fig. S3(b)), about 3 times 
smaller than our AIMD value of (9.6±2.0) x 10'^ Nsm“^. 
This is because the optB88-vdW functional predicts a 
larger corrugation of the free energy compared to the 
force field used |2S](see Fig. S4). Since optB88-vdW 
overestimates the absolute adsorption energy of a water 
monomer on graphene when compared to benchmark dif¬ 
fusion Monte Carlo and Random Phase Approximation 
results HH]) we cannot expect optB88-vdW to capture 
the corrugation of the free energy with absolute preci¬ 
sion. However, as explained earlier, because this func¬ 
tional successfully captures various other properties of 
water, graphite, BN and water at interfaces we expect it 
to predict the correct relative water/graphene and wa- 
ter/BN interaction strengths and hence the correct in¬ 
crease in the friction on BN. Since the increase in fric¬ 


tion is not captured by force field molecular dynamics 
(see Fig. S3(b)), the findings here stress the importance 
of accounting for electronic structure effects when inves¬ 
tigating transport properties at complex liquid/solid in¬ 
terfaces. Nevertheless, it is possible that improving the 
description of the force fields to include polarizable mod¬ 
els and partial charges on B and N may reproduce the 
observed friction increase. For instance, it has been found 
that including polarization effects in force field molecular 
dynamics studies has an effect on the diffusion of liquid 
water on charged graphene [34j . 

Finally, there has been increasing interest in connect¬ 
ing wetting properties to the friction coefficient (see e.g. 
Ref. [38]) and a relation between the slip length and the 
contact angle has been found to hold for a wide num¬ 
ber of liquid/solid interfaces. Here we have seen that the 
structure of liquid water on graphene and on BN is strik¬ 
ingly similar and previous ab initio work reported also 
contact angles of 86° and 87° on graphene and BN, re¬ 
spectively m- Yet, friction is about three times larger on 
BN, highlighting that a simple dependence on the wetting 
properties cannot be established for these two systems. 
Instead, we have demonstrated a dependence of the fric¬ 
tion on the free energy of the water contact layer. In sys¬ 
tems like graphene and BN, where the corrugation of the 
potential felt by the water is not directly proportional 
to the water/solid interaction strength, the free energy 
profile provides a closer estimate of the potential energy 
landscape corrugation. Possible other examples are car¬ 
bon nanotubes, which like graphene, have been shown to 
depart from the scaling law that relates the friction to the 
water contact angle m-, and most likely BN nanotubes 
and other van der Waals layered materials m- With the 
aim of designing nanofluidic devices which exhibit fric¬ 
tionless fluid transport, materials may be engineered to 
favor a smooth potential energy landscape independently 
of the fluid/solid interaction strength. In this manner for 
instance, the friction may be tailored whilst maintaining 
the same wetting properties. We have shown that be¬ 
cause of electronic structure effects, friction is larger on 
BN. Limiting the magnitude of such interactions between 
the liquid and the substrate using for instance homopolar 
surfaces and avoiding the presence of H-bonds is benefi¬ 
cial for the design of smooth interfaces. Further, com¬ 
pressive or tensile stress can be applied to a substrate to 
favor a smooth contact layer [48] . 

In conclusion, we have reported on extensive ab initio 
molecular dynamics studies of liquid water on graphene 
and BN. In so doing we have tried to bridge the gap be¬ 
tween the molecular structure and energetics of water on 
layered materials and the complex transport of water at 
the nanoscale fully from ab initio methods. A comparison 
between the two systems reveals that while the structure 
of the liquid film is very similar, slight differences in the 
water contact layer not related to the wetting proper¬ 
ties of the interfaces give rise to a remarkably different 
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water slippage on the two sheets. The three-fold fric¬ 
tion increase on BN is induced by a larger corrugation of 
the energy landscape compared to graphene, because of 
more pronounced electrostatic and exchange interactions. 
Overall, this work paves the way for the study of trans¬ 
port properties in yet more complex liquid water/solid 
interfaces using ab initio methods. For example, sys¬ 
tems where water is liable to dissociate or where ions, 
defects, or external electric fields are present could all be 
examined. We also hope that our work will stimulate the 
study of nanoscale water friction on BN-nanostructures 
and other layered materials using e.g. atomic force mi¬ 
croscopy [39]. 


SUPPLEMENTARY INFORMATION 

Further computational details and tests on the struc¬ 
ture, energetics and friction coefhcient of the wa¬ 
ter/graphene and water/BN interfaces. 
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